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Abstract 



The scaling properties of waves of topplings of the sandpile model on a Sierpinski gasket are 
investigated. The exponent describing the asymptotics of the distribution of last waves in an 
avalanche is found. Predictions for scaling exponents in the forward and backward conditional 
probabilities for two consecutive waves are given. All predictions were examined by numerical 
simulations and were found to be in a reasonable agreement with obtained data. 



1 Introduction 



Sandpile models form the paradigmatic examples of the concept of self organised criticality (SOC) 
P], ||]. This is the phenomenon in which a slowly driven systems with many degrees of freedom 
evolves spontaneously into a critical state, characterised by long range correlations in space and 
time. 



In the past decade much progress has been made in the theoretical understanding of sandpile 
models. This is especially true for the Bak-Tang-Wiesenfeld (BTW) model [Q, g], where, following 
the original work of Dhar ||, a mathematical formalism was developped [4-8] that allows an exact 
calculation of several properties of the model such as height probabilities ||, |5[ , the upper critical 
dimension [|| and so on. 



Despite all this work, it has however not been possible yet to give a full and exact characterisation 
of the scaling properties of the avalanches in the BTW-model. In recent years it has become 
increasingly clear that, especially in two dimensions, avalanches are to be described by a full 



multifractal set of scaling exponents |9|, 10]. This spectrum of exponents has been calculated with 



high numerical precision, but at this moment there is no clue how it can be determined by an 
analytical approach. 



Avalanches can be decomposed into simpler objects called waves [Q]. The probability distribution of 
waves seems to obey simple scaling [ pd| ] and the exponent describing that scaling is known exactly, 
both for the general wave |llj and for the last wave of each avalanche ||]. Since in dimensions 
d > d c = 4, multiple topplings are extremely rare, it is to be expected that in these situations, 
wave and avalanche statistics obey the same scaling properties. 



More recently, the distribution of two consecutive waves has received considerable attention [12, 13 1. 
The conditional probability that the k + 1-th wave has size Sk+i given that the previous wave had 
size Sfc, P(sfc_j_i|sfc), is the first quantity to study when one is interested in correlation effects in 
waves. It are these correlations that make waves and avalanches different. Paczuski and Boettcher 



] 12| proposed, on the basis of extensive simulations, that P(sk+i\sk) has a scaling form 



P(s k+ i\s k ) ~ s^F(^±k) (1) 

where for large x, F(x) ~ x~ r , while F(x) — > constant for x —* 0. Numerical estimates for these 
exponents in d = 2 are (3 ~ 3/4, r ~ 1/2. At this moment, no exact values for these exponents 
are known. In a very recent work, Hu et al. ]l4| study the 'backward' conditional probability 
P(sk\sk+i) and find that it obeys a similar scaling law 

P(s k \s k+1 ) ~ sfF(^) (2) 

Sk+l 

where for large x, F(x) ~ x~ r , while F(x) — > constant for x — > 0. These authors give arguments 
that show that 

P + r = T lw (3) 

where T\ w is the scaling exponent describing the size distribution of the last wave. The relation 
(01) is consistent with the numerical data for the square lattice where both Euclidean and fractal 



2 



dimensions of waves are 2. At the same time, it is desirable to get an independent verification of 
(|3|) using lattices of different dimensions. 



In the present paper we study the properties of waves on the Sierpinski gasket, continuing previous 
work [15|. Using the methods of analysis introduced in 0, we obtain precise numerical estimates 
for the exponents Ti w ,(3,r,(3 and r. Also in this case equation (||) seems to be well satisfied which 
indicates that short time correlations in waves admit an analytical treatment. 



2 The sandpile model on a Sierpinski gasket 



The BTW sandpile model can be defined on any graph, but for definiteness we will introduce it 
in the context of the Sierpinski gasket (see figure 1). Each vertex (apart from the three boundary 
sites) of this graph has four nearest neighbours. To each such vertex i we associate a height variable 
Zi which can take on any positive integer number. We also introduce a critical height z c , which we 
will take equal to four for all vertices. The number of sites in the lattice, N, is trivially related to 
the number of iterations n used in constructing the fractal. The dynamics is defined as follows. On 
a very slow time scale we drop sand at a randomly selected site and thereby increase the height 
variable by one: z% — > z% + 1. When at a given site, z% > z c , that site becomes unstable and topples 

z j ~^ z j ~ 

where 

4 i = j 

Aj = —1 if i and j are nearest neighbours (4) 
otherwise 

Through toppling, neighbouring sites can become unstable, topple themselves, create new unstable 
sites, and so on. This avalanche of topplings proceeds on a very fast time scale and no new grains 
of sand are added before the avalanche is over. Sand can leave the system when a boundary site 
topples. An avalanche is over when all sites are stable again. 

For further reference it is also necessary to introduce the matrix G, called the lattice Green function, 
which is the inverse of A. 

It is not difficult to see that the order in which unstable sites are toppled does not influence the 
stable configuration which is obtained when the avalanche is over. This Abelian nature of the 
sandpile model allows the introduction of the concept of waves |jj which are introduced in the 
following way. Suppose that an avalanche starts at a site i$ and that after a few topplings iq 
becomes unstable again. One can then forbid iq to topple again and continue with the toppling of 
other unstable sites untill all of them are stable. It is easy to show that in such a sequence all sites 
topple at most once. This set of topplings is called the first wave. Next, we topple the site iq for 
the second time. If after some topplings it becomes unstable again, we keep it fixed, and topple 
all the other unstable sites. This set of topplings constitutes the second wave. We continue in this 
way untill we finally reach a stable configuration. We can in this way decompose any avalanche 
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Figure 1: Sierpinski gasket with n = 2. The arrows indicate the three sites where sand can leave 
the system. 



in a set of waves. The probability distribution P w (s) that an arbitrary wave involves s topplings 



obeys simple scaling [11] 



P w (s) « (5) 
(For the moment we neglect finite size effects which will be taken into account in section 4.) 

Because in a wave, sites topple at most once, waves are simpler objects to analyse than avalanches. 
Without going into details we summarize the following important properties of waves [|J 

• There is a one-to-one correspondance between waves and the two-rooted spanning trees on a 
graph which consists of the Sierpinski gasket and one extra site called the sink. The sink is 
connected with two edges to each of the three boundary sites of the Sierpinski gasket. 

• The element Gij of the Green function is given by the ratio of the number of two-rooted 
spanning trees (in which i and j are in the same subtree) to the number of one-rooted 
spanning trees. Moreover Gij is also equal to the expected number of topplings at site j when 
a grain of sand was dropped at site i, which is proportional to the probability that a wave 
started at % reaches j. 



From these results, it follows that P W (R) = dG{R) / dR where R is the linear size of the wave. The 
asymptotic behaviour of the Green function on an arbitrary lattice is 

G{R) rsj R d ^- d f 

where df is the fractal dimension of the lattice, and d w the dimension of a random walk on the 
lattice. Therefore, P W (R) ~ R">w—df-i _ p rom definition of fractal dimension, s ~ R f we finally 
obtain 

Pw(s) ~ S d ™/ d !~ 2 (6) 
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so that 



r w = 2 - <± (7) 
d f 

a result first derived in |i~5[ j. For the particular case of the Sierpinski gasket, d/ = log 3/ log 2, = 
log 5/ log 2, so that t w = log (9/5) / log 3 » 0.535, a result which is nicely consisted with the available 
numerical data. 

The properties of the last wave in a given avalanche are of special interest for us (see section 3). 
The probability distribution Pi w (si w ) that the last wave has si w topplings obeys the scaling law 



HI 



PUs lw ) ~ S/ ' Tiw (8) 

From the definition of waves it follows that the last wave has the property that the site iq is on the 
boundary of the wave. Let us denote by ds the fractal dimension of the boundary of an arbitrary 
wave. The number of points on the boundary of a wave of size s is then of order s dB ^ df . Hence, 
the probability that a given site is on the boundary of a wave of size s is proportional to s c( s/ c( / _1 ) 
which should also be proportional to Pi w (si w )/P w (s). Using (6) we immediately obtain 

d w + ds , , 

n w = 3 (9) 

d f 

In two dimensions, d w = 2 and ds = z = 5/4 where z is the fractal dimension of the chemical path 
on a spanning tree. One thus finds ti w = 11/8 ||. On a fractal the relation between ds and z may 
be more complicated. In fact, it was shown in [|16| that for a deterministic fractal 

ds = z - d w + df (10) 

so that, from (^) we obtain 

T~lw = 2 -j— (11) 
d f 



The exponent z on the Sierpinski gasket was also calculated in [16] with the result 



z = log [(20 + v205)/15]/ log 2. Thus one finally obtains for the case of the Sierpinski gasket 

log[(20 + V205)/135 U l247 
log 3 

In section 4, we will present numerical estimates for t\ w that are fully consistent with this prediction 



3 Distribution of consecutive waves 



In order to characterise the statistical properties of waves more completely, it is necessary to go 
beyond the description on the basis of the distribution P w (s) only. Following the work of Paczuski 

lr The derivation of ri w in |l5[ used implicitly that the graph is selfdual, which is correct for the square lattice but 
not for the Sierpinski gasket. This error was pointed out by one of us (VBP). The correct result is that given in (|12[). 
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and Boettcher 12] we now turn to a study of the conditional probability P(s k +i\s k ) where s k is 
the size of the fc-th wave. If the size of consecutive waves is a Markov process, this conditional 
probability is sufficient to describe the evolution of wave sizes. In figure 2. a, we show our data 
for this conditional probability for Sierpinski gaskets with n = 9 (N = 29526). All our data were 
obtained by studying at least 1000 x N avalanches. The figure shows the best fit of our data to the 
scaling form (|l|) proposed by Paczuski and Boettcher. Unfortunately, it is not possible to obtain 
very accurate estimates for the exponents (3 and r in this way. We will come back to this issue in 
the next section. 



Recently, it was pointed out that the 'backward' conditional probability P{s k \s k +i) is also of interest 
because it is possible to relate the exponents (3 and f (see (g)) to t\ w . In figure 2.b we present our 
data for this quantity, again for the case n = 9. 



We now repeat briefly the argument given in []14|]. We begin by rewriting (||) in a normalised form 

%l%i)~(f)"^(f) «j& (13) 



Let's next consider the situation in which s k 3> s k+ i so that the argument of F in ( |l~3| ) is large. In 
that case it must be so that the (fc + l)-th wave has a non-empty intersection with the boundary of 
the previous wave (see figure 3. a). Indeed, assume the opposite so that the (k + l)-th wave covers 
a small region inside the much bigger fc-th wave (figure 3.b). But such a situation is forbidden, 
since all sites inside the k-th. wave return to their original height after the wave has passed. The 
(k + l)-th wave must therefore follow the motion of the previous wave untill it hits the boundary 
of the fc-th wave from where it can follow a different evolution. Therefore, the situation of figure 
3.b cannot occur. 



Then, consider figure 3. a on a coarse grained scale by performing a rescaling of the order of Rk+i, 
the linear size of the k + 1-th wave. On that scale, the geometry of figure 3. a resembles that of last 
waves with s k+ i playing the role of the origin of the avalanche, and s k that of the last wave. Hence 
we arrive at the conclusion that in the limit s k ^> s k+ i, the distribution of s k coincides with that 
of the last wave. From (|l^), @ and the asymptotic behaviour of F, the equality (0) then follows. 

Once this result has been obtained it is possible to obtain also a relation for the exponents (3 and r 
that appear in the scaling form ([l]). The joint distribution P(s k , Sfe+i) can be written in two ways 
using either the forward or backward conditional probability 

P(s k ,s k+1 ) = P(s k \s k+1 )P w (s k+1 ) 
= P(s k+1 \s k )P w (s k ) 

We then insert (13), © and a properly normalised version of ([!]) and get 

In the case s k S> s^+i we insert the proper limiting behaviours of the functions F and F, and 
immediately obtain, using (0) 

(3 = 1 + T w - T iw (15) 
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Figure 2: Scaled conditional probabilities P(sk+i\sk) (a) and P(sk\sk+i) (b). The results were 
obtained for a Sierpinski gasket with n = 9. 
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a 



b 



Figure 3: When Sk+i <C Sfc the k + 1 the wave must have an intersection with the boundary of the 
k-th wave as shown in a. The situation in b cannot occur. 

A final equality between exponents can be obtained by investigating the limit st -C Sfc+i in ([HI). 
Inserting the appropriate scaling behaviours one obtains 

P + r =l + Tw -P (16) 

In d = 2, ( |i~5| ) leads to the predictions j3 = 5/8 and /3 + f = 11/8. The value of f3 is not too far 
from the numerical estimate (3 ~ 3/4 reported by jl2|1 ) while in |l4| numerical evidence is presented 
that is in agreement with the prediction (^). In the following section, we investigate the situation 
on the Sierpinski gasket. 



4 Numerical results 

In order to analyse our data we have used the method introduced in || , in which one investigates 
the moments {s q ) of the distribution P(s,L), where we now explicitly take into account the size L 
of the system. In our case, L = 2 n . If P(s,L) has a simple scaling form 

P(s,L) ~ s- T H(s/L d f) (17) 

these moments should be proportional to simple powers of L, (s q ) ~ L^") wri ere a(q) = d f (l-T + q) 
for q > t — 1, and a(q) = for q < r — 1. 

In principle, an analysis of the moments is most instructive when one is interested in the presence 
of possible multifractal scaling (instead of simple scaling). In that case, the function a(q) becomes 
nonlinear. However, even in the absence of multifractality, this kind of analysis has many advan- 
tages. In the case of Sierpinski gasket, the probability distributions for the size of avalanches, waves 
and last waves show strong oscillations superposed on the pure power laws (see the figures 2,4 and 
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Figure 4: The moment exponent cri w (q) for last waves. The straight line gives a best fit to the high 
(/-data, from which t\ w can be obtained. 



5 in pq| ) . This is a consequence of the discrete scale invariance of the system. These oscillations 
make the determination of the scaling exponents a hard task. The moments {s q ) have the big ad- 
vantage that they are averages over the distribution and hence the effects of the oscillations almost 



completely disappear. If one then assumes simple scaling, as we know is correct for waves [11|, one 
can further reduce any remaining fluctuations by fitting the values of a(q) (for q big enough) to a 
straight line. The slope of the line should equal df, and the intersection with the g-axis gives r — 1. 

We tested this method of analysis for the last wave. We performed extensive simulations for 
Sierpinski gaskets with 5 < n < 9. From the statistics of the last waves we could then estimate 
cri w (q). The results are shown in figure 4. The curvature at low q is a finite size effect. From a fit 
of our data in the regime q > 1.5 we obtain the estimate T[ w = 1.28 ± .03. This is very close to the 
exact result which we obtained in section 2 0. 

We have followed a similar scheme of analysis for all the other exponents introduced in section 3. 
Take as a concrete example the sum of exponents (3 + f. This can be obtained as follows. For 
Sk 3> Sfe+i, P + {sk) = P{sk\sk+i) ~ s fc r ~^- Instead of analysing this distribution itself we look at 
the moments of P + (sf c ) which are expected to scale as L^^. A plot of (J + {q) is shown in figure 
5. From an analysis of the data in this figure we obtain the estimate (3 + f = 1.24 ± .01, which is 
very clearly consistent with the prediction (||). 

By analysing the small Sk behaviour of P{sk\sk+i) in a similar way, we can estimate (5 = .16 ±0.05. 

2 This value of ti w is also much lower then that obtained in jlB) from a direct fit to (|^). We are currently performing 
a multifractal analysis of data for the size and the area of avalanches in the BTW (and in a stochastic sandpile) 
model, on the Sierpinski gasket. The results will be published elsewhere. 
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Figure 5: The moment exponent cr + (q) (see text). The straight line gives a best fit to the high 
g-data, from which /3 + f can be obtained. 

Continuing in this way for the forward conditional probability P(sk+i\sk), we find from the large 
Sfc+i-behaviour (3 + r = 1.30 ± 0.01, while from the data for small Sfc+i we finally obtain (3 = 
0.35 ± 0.05 (see figure 6). This numerical result is not too different from the prediction following 



from (15) which gives (3 = .288. 



Finally note that also the relation ( |16[) is rather well satisfied 



5 Conclusions 



In this paper we investigated the properties of waves in the sandpile model on a Sierpinski gasket. 
We gave predictions for the exponent describing the last wave in an avalanche and for the scaling 
exponents occuring in forward and backward conditional probabilities for consecutive waves. These 
predictions were tested by extensive simulations and were found to be in good agreement with the 
numerics. 

Results such as those shown in figure 5 and figure 6 also show no clear evidence for any multifrac- 
tality which would show up as a curvature in the plots of a(q) for big enough q ([0 ). 

We are currently investigating the presence of multifractality of avalanches on the Sierpinski gasket. 
If, as is the case in d = 1 ( [|17]] ) and in d = 2 ( ||, [l(| ) , such multifractality shows up , the interesting 
question arises how such a phenomenon can be built up on the avalanche level when it is absent at 
the level of waves and consecutive waves. 
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Figure 6: The exponents of the moments of P(sf t +i\sf (: ) for s^+i -C s^. The straight line gives a 
best fit to the high g-data from whicht the exponent (3 can be estimated 
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